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Abstract 



Coupled PDE models describing both morphological and chemical evolution, such as recently 
introduced by Bradley and Shipman, offer a promising potential explanation for the highly-ordered 
patterns sometimes observed during the irradiation of binary compounds with low-energy ions. 
However, as currently employed in several contexts, these models require specific conditions on 
a number of unknown parameters to explain this behavior. Specifically, they require that the ef- 
fect on surface morphology of atoms sputtered away from the surface is larger than the effect of 
the atoms merely redistributed to new locations on the surface, a formerly-common assumption 
that has recently been overturned for several pure materials. Here, we develop a framework in 
which many of the parameters of the BS theory can be estimated using the results of molecular 
dynamics simulations. By adapting the recent theory of "crater functions" to the case of binary ma- 
terials, we are able to establish relationships between the parameters in coupled PDE models and 
various moments of the statistically-averaged crater function associated with single ion impacts. 
Applying our framework to the commonly-studied GaSb system, we find that the requirement of 
erosive dominance is not fulfilled, motivating the consideration of other potential pattern-forming 
mechanisms. 
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I. INTRODUCTION 



The irradiation of flat surfaces by low-energy ions can lead to the emergence of ordered 
arrays of nanoscale surface features IfTH , and offers a potential to assist in the fabrica- 
tion of nanostructured devices with regular structure U2H . Initially observed on a wide 
variety of target materials Q3H5B , careful experimental investigation has led to a growing 
consensus that these structures only appear during the irradiation of binary materials, or 
environments in which one or more materials is simultaneously co-deposited during the ir- 
radiation process ll6l- [TTll . This consensus has motivated a number of theoretical treatments 
of irradiated binary materials, tracking the evolution of both the morphology and the con- 
centration. Extending the pioneering work of Shenoy Chan, and Chason IH2H . Bradley 
and Shipman (BS) have recently introduced such a theory exhibiting the first physically- 
grounded explanation for the strong ordering of the observed nanostructures ITT31 |14ll . 
This result has sparked great interest, inspiring additional studies on normal-incidence ir- 
radiation of two-component materials lH"5H . including a variant where the second material 
is not originally present, but instead co-deposited simultaneously with the ion irradiation 
|fl"6ll , and a similar framework for the study of ion-assisted deposition, where both materi- 
als are deposited along with the ions, and the deposition rate exceeds the sputtering rate 

eh. 

Despite its success, however, the family of models initiated by Bradley and Shipman 
face some important open questions. At a general level, they contains many parameters 
whose values are so far undetermined. Because the BS theory only admits ordered struc- 
tures for specific regions of parameter space, it is important to identify whether common 
experimental systems have parameters within those regions. So far, such estimation does 
not appear to have been performed. More specifically, it assumes that the fundamental 
driving force behind the instability is the erosive mechanism identified by Bradley and 
Harper [18], implying that the surface should be unstable even in the absence of a second 
species. This has formerly been a common assumption in the field, but has recently been 
overturned for the case of the pure materials Silicon and Germanium llT9ll20ll . where flat 
surfaces are observed at normal incidence, attributed to the dominance of redistributed 
atoms over sputtered atoms [|2T1 - [23l (though beam-induced stress within a thin layer of 
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amorphized material also serves to stabilize surfaces at normal incidence ||24l - l26l0 . The 
important question we will pursue, then, is to determine whether this dominance extends 
to the case of binary materials, and if so, to identify the consequences for existing theories. 

We will begin in Section II with a necessary, brief review of the key features of the 
Bradley-Shipman theory: its equations, parameters, the ordered patterns it predicts, and 
especially the parameter regime necessary to obtain those patterns. Then, in Section III, 
we begin our main task: the development of a framework for the estimation of several 
of the parameters, by means of an extension of the theory of crater functions [|23l 12711 
to the case of binary materials. We demonstrate how this generalization leads, as in the 
case of pure materials, to a continuum equation for the evolution of the surface height in 
terms of moments of the crater functions. Comparing that equation to those of Bradley 
and Shipman, we then obtain relationships between parameters in the latter theory to the 
crater function moments in the former, which can be estimated using molecular dynamics 
simulation. Finally, in Section IV, we perform such simulations of irradiated Gallium Anti- 
monide (GaSb) - the experimental system in which ordered dots were first observed Ifflh 
and which remains of experimental interest within the ion-beam patterning community 
Il28l |29ll . Using the the framework developed in Section III, we estimate parameters for 
this system, finding that redistributed atoms continue to dominate over sputtered atoms 
for GaSb under low-energy argon irradiation, and concluding that the BS model may not 
actually be able to explain the finite wavelength bifurcation in this system. A discussion of 
the implications follows in Section V. 

II. REVIEW OF THE BRADLEY-SHIPMAN THEORY: MODEL AND STABILITY ANALYSIS 

We begin by briefly summarizing the recent model of Bradley and Shipman UT3l IT4H 
for irradiated binary materials, which itself builds upon the earlier model of Shenoy et al. 
IIT21I . Because extensive analysis of this class of models exists elsewhere, both in general 
H30H . and in the context of ion irradiation [fl2l - [T4l IT7H , we will be as brief as possible. 
However, because our aim is to test specific aspects of this theory and others like it, we must 
include those results that are directly relevant to the molecular dynamics simulations we 
perform at the end. (Note: when referring to specific equations within the BS theory, we 
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refer to the longer and more detailed Ref. [14], despite the earlier publication of Ref. [fT3tl ~). 



A. Model 

The Bradley-Shipman model tracks the evolution of a height field z = h(x, y, t) de- 
scribing the irradiated surface, and concentration fields ca {x,y,t) and c B (x,y,t) of two 
components A and B. Under the effects of preferential erosion, a steady state is reached 
in which the material is receding with with constant velocity v and constant surface con- 
centrations CAfl and cb ;0 of A and B atoms, respectively. The stability of this steady state 
is investigated by introducing perturbations to this steady state via the equations 

h = —v t + u (x, y, t) 
ca = c A fl + <f)(x,y,t) 
c B = c B ,o + (1 -<j>{x,y,t)) 

where u {x, y, t) describes the perturbation to the height field, and <fi (x, y, t) describes the 
perturbation to the concentration field of species A. After some significant modeling and 
analysis (Eqs. (3)-(13) of Ref. [14]) the following linearized equations for the evolution of 
u and (f) are obtained: 



(/II 

^ = -A<f> + CV 2 u - DV% (1) 

at 

^ = -A'<f> + BV 2 <P + CV 2 u , (2) 

at 



where we have continued to employ the reasoning in Ref. |fT4ll used to neglect terms BV 2 c 
in Eqn.Q and 5V 4 h in Eqn.(|2]). In a slight deviation from Ref. []T4l , we will assume that 
species A is preferentially sputtered (rather than species B), and define coefficients in 
such a way that zeroth- and fourth-order terms in Eqs.([T])-(|2]) have negative sign, whereas 
second-order terms have positive sign. The coefficients in Eqn.Q for the height evolution 
are defined via 
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A = P n [A' A (c Afi ) - A' B (c B ,o)] > 

C = Q [(p A (c A ,o) + Hb (cb,o)) - a (A A (c A ,o) + A B (cb,o))) (3) 

D = [c A , ^ + c B>0 D B ] -f-^=? > 
and in Eqn.Q for the concentration evolution, via 

A' = [c B ,bA' A (ca,o) + c A)b A' B (c B ,o)} > 
n Q 

B' = ^~ [c B , b D A + c Ab D B ] > (4) 

C = ^ [cb,6/44 (ca,o) - Ca,6A*B (cb,o)] 

Here P is the power deposited by the ions per unit surface area on a flat surface, is the 
atomic volume (taken to be the same for both species), n s is the total number of mobile 
surface atoms per unit area on the surface, 7,, is the surface energy, k B is Boltzmann's 
constant, T is the temperature, A is the amorphous film thickness, c Ai0 and cb,o are the 
concentration of A and B atoms in the film at steady state, c Atb and cs,b are the concen- 
tration of A and B atoms in the bulk, D A and D B are the diffusivities of A and B atoms, 
(c A ) and A B (c B ) are proportionality constants linking the deposited power to the sput- 
tering rate, and /i A (c A ) and jj B (cb) are proportionality constants describing preferential 
redistribution of material, an effect first proposed in Refs . 11131 [T4l and playing a critical 
role therein. 



B. Generic Stability Analysis 

To determine the presence or absence of an instability, is is sufficient to consider a one- 
dimensional perturbation, because of the isotropic nature of normal-incidence irradiation. 
Without loss of generality, we may orient this instability in the x-coordinate direction, 
giving 



u (x,t) 




«1 


ikx+at 


(5) 


OM) 
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where k is the wavenumber of the perturbation, a the growth rate of that perturbation, and 
ui and 0i are undetermined constants describing the relative phases and magnitudes of the 
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height and concentration modulations. Inserting the perturbation ([5]) into the governing 
equations ([T])-(|2]), we obtain, in matrix form, 

a + Ck 2 + Dk A A 

C'k 2 a + A' + B'k 2 

For a solution [ui, (f>i\ of this equation, the determinant of the matrix must be zero, giving 
a the dispersion relation a (k) in the quadratic form 

a 2 + to + A = 0, (7) 

where r (A;) and A (k) are respectively the trace and determinant of ([6]) when a = 0: 

r (k) = A' + (C + 5') A; 2 + Dfc 4 (8) 
A (k) = {CA 1 - AC) k 2 + (CB 1 + DA') k 4 + (DB') k 6 '. (9) 

For the purpose of identifying stability properties, it is sufficient to consider only the root 
associated with the '+' sign in the quadratic formula which, following common practice, is 
denoted a + (k): 

a + (k) = \(-r + vV - 4A) . (10) 
This gives the growth rate of the faster-growing solution when the discriminant is positive. 

C. Discussion: Requirements for Ordered Patterns 

The expected behavior of an initially-flat surface is determined by nature of the func- 
tion er + (k), extensively discussed in Refs. fl30l I3TT1 . Wavenumbers k for which a + (k) < 
are stable and decay over time, while those k for which a + (k) > are unstable and grow 
over time. The essential first ingredient for producing ordered patterns is the presence of 
a narrow band of unstable wavenumbers surrounded on both sides by stable regions. In 
particular, the very-smallest wavenumbers corresponding to very long wavelengths must 
be stable; otherwise the system will tend to exhibit disordered kinetic roughening phenom- 
ena. An important accomplishment of the Bradley-Shipman theory was its identification of 
such a narrow band in the Equations ([l])- ([2]). Although Eqs.Q-Q contain six parameters, 
the conditions under which such an instability can occur are quite limited. This is most 
easily seen via a brief qualitative simplification of the analysis in Ref. Ifl4ll . 
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01 



0. 



(6) 



First, we consider the conditions under which an instability of any kind may exist. 



Examining Equation ( [10] ), we see that the system is stable if both r > and a > 0, but 
unstable if either one is negative. From the definitions @-((9]), we see that because the 
parameters A, D, A', and B' are all positive, this can only occur if one of either C or 
[A'C — C'A) is negative. Second, turning to the requirement, for ordered patterns, that 



long wavelengths be stable (to form a narrow band), we examine Eqs. (|8|)-(|10j) in the limit 
as k — > 0, and find that for small values of k 

lim a+ (k) ph {A'C - C'A) k 2 + O (k 4 ) . (11) 

To stabilize the longest-wavelength modes, we must therefore impose the requirement 



A'C > C'A. (12) 

From the above considerations, we immediately conclude that an instability with a nar- 
row band of unstable wavenumbers is only possible if (a) C is sufficiently negative to drive 
the instability, and (b) 7 is sufficiently negative to stabilize the long waves. These are strik- 
ing requirements of the Bradley-Shipman instability. The first implies the existence of an 
morphologically-driven instability IfTBl even in the absence of concentration modulations 
(ci = 0) . While formerly a common assumption, this has has recently been overturned for 
pure Silicon and Germanium |fT9l |20~H . where there is strong evidence that in fact C > 
due to the dominance of redistributive over erosive effects [23]. Although binary systems 
are chemically different from pure silicon, the collision cascade leading to both erosion 
and redistribution is a predominantly kinetic process (although chemical effects have been 
identified in collision cascades, they do not affect atom redistribution strongly H32l 13310 . 
Hence, the dominance of redistribution over erosion might be expected to persist in bi- 
nary system. The second requirement of the BS instability, of strongly negative C , is also 
somewhat surprising, for it implies that the preferentially-sputtered material is much less 
likely to be displaced within the collision cascade. These features highlight the importance 
of being able estimate parameters within the BS theory, and we now illustrate how this 
may be done via molecular dynamics. 
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III. COEFFICIENT ESTIMATION VIA MOLECULAR DYNAMICS 



A. Crater Functions for compound targets 

To estimate parameters in two-component models such as the BS theory, we must first 
develop a framework within which the results of MD simulation may be connected to the 
coefficients of Eqs.([T])-([2]). To this end, we will generalize the "crater function" theory 
for surface morphology evolution Il23l |27ll to the case of a multi-component material, 
obtaining continuum equations of a form similar to Eqn.Q. Essentially, we need only 
revise the definition of the "crater function" Ah(x,y; 9), describing the average change in 
surface morphology to a flat surface due to the impact of a single ion at the origin, so that 
it is composed of four parts: 

Ah (x, y- 6, c A , c B ) = Ah c ™ s - + Ah e ™ s - + Ah v f lst - + Ah r f lst - (13) 

This merely states that the change in surface topography due to a single ion impact is 
the superposition of effects due to both eroded and redistributed atoms, of both Ga and 
Sb. A central benefit of this approach is that any subsequent analytical treatment of Ah 
treats all these effects identically, and so long as all operations on Ah are linear, the final 
contribution to continuum equations can be broken out into component parts at the end 
of the analysis. Hence, because the analysis of Ref.[27] (and a simplified approximation 
for flat targets in Ref.[23]) has already been performed for a general Ah, all the results 
therein continue to be valid for the present case. 

Using Eqs.(46) and (48) in Ref. U27H as a template, we define 

AhT s - = -n A £ a (xf) 

i 

Ah%° a - = -n B J2H^) 

j 

Ahf lst - = n A ^5 (xf ) - 5 (xj) 

k 

Ah^-= fi B ^r<5(xO-5( x 

I 

where the denote atomic volumes, i,j denote sputtered A,B atoms, k, I denote displaced 
A,B atoms, and and x F denote initial and final target atom positions associated with 
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an ion impact. When the crater function is integrated to obtain moments, these four 
components produce associated contributions to the zeroth and first moments. Following 
Eqs.(47) and (49) of Ref.H27t we obtain: 



hlf 



M 



(0) 
B 



M 



(1) 

A,eros. 



M 



(i) 



M 



B,eros. 



(1) 

Ajredist. 



M 



(i) 

S.redist. 



Ah^-dA 



Ah c ™ s -dA 



xA/i r B edist ^ 



i 

j 

i 

i 

^(xf-xj) 



(14) 



Finally, following a multi-scale analysis, these contributions to the moments produce asso- 
ciated terms in the governing equation for the height field 



ht « (Y A + Y B ) 



+ [ S X ' + b y + by + by' ] Al, 



X 



(15) 



+ [Sy' + Sy' + by + by' ) k 



yti 



with definitions given by Eqn.(4) of Ref. H23ll 



y z {9, a) 



I Q cos (0) M^ 0) (0, c, 



Q Z,type 



d_ 
d6 



sp type {e, a) 



I cos(6)M^ type (e,c t ) 



'l co S (9)cot(e)M^ ype (e,c i ) 



(16) 



where Z is either A or B and 'type' is either 'eros.' or 'redist.' (see Section [IV C| for further 
comments on these formulas). In addition, J is the flux through a plane perpendicular to 
the beam, which is not defined within the simulation of a single impact. 



B. Relationship to BS Theory 

Upon comparison with BS theory, it is easily seen that the coefficients of the terms 
h xx and h yy (which are equal at normal incidence due to rotational symmetry) correspond 
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exactly with the coefficient C in Eqn. ([!]). However, the various quantities Y and S in Eq|l5"1 
allow estimation of the physical parameters underlying C, which critically also appear in 
other parameters. We will see that the results of molecular dynamics in the form of the 



coefficients 16 provide enough information to estimate fully four coefficients {A,C, 0,7} 
from the Bradley-Shipman theory. 

Before we proceed, we make one important simplifying assumption. The various com- 
ponents of the moments and are functions not only of angle, but also of the 
concentrations of the various species within the target. Similarly, within the BS theory, 
the quantity Aj - associated with sputter yield - and presumably also ^ - describing the 
redistributive flux - share this dependence. [The latter dependence was not discussed in 
Ref.[14], but presumably exists.] As an alternative to performing simulations at many 
different concentrations, we will here make the assumption on the moments that 

<' (cs P .) = M<°> (I) x (2 Ci ) 

)l\ (17) 

M« (c sp .) = AfW I -J x (2c,) 

and the corresponding assumption, within the BS theory, that 

A,- (a) = Ai I - I x 2 Ci 



A (18) 

That is, each species both sputters and redistributes at a rate linearly proportional to its 
concentration in the film, with rate zero when its concentration is zero, and rate equal 
to our measured rate when its concentration is one half. This will allow estimation of 
derivatives with respect to concentration in the vicinity of the naturally-occurring 50/50 
concentration, without excessive numerical simulation (see Section [TV C| below for further 
discussion) . 

With the approximations given by Eqs.([l7])-(18), we are now able to obtain a direct 



relationship with the theory of Bradley and Shipman by directly comparing Eqn.(fl5|) with 
Eq.(7) in Ref. ITT411 ■ Although these equations were obtained by different means, they de- 
scribe precisely the same phenomenon, based on an identical set of physical mechanisms. 
We therefore equate them in a term-by-term manner for the case of normal incidence 
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(6 = 0). Using also Eqs.(3-6) of Ref. [fT4ll , we see that the various coefficients ( [16] ) ex- 
tracted from the crater function moments must be equal to parameter combinations in 
Ref. [14] as follows: 



(19) 



n A P A Afl = -2Y A (0,0.5) x c A)0 
n B P A Bfi = -2Y B (0,0.5) xc Bfi 
n A aA Afi = -2S^ os - (0,0.5) x c Afi 
n B aA Bfl = -2S x ^ os - (0, 0.5) x c Bfl 
^A^Afi = 2Sx'™ dmt - (0, 0.5) x c Afi 
^Bf^B,o = 2SxY dmt ' (0, 0.5) x c Bi0 



where the Y and S terms are estimated for 9 = and at the 50/50 concentration, and the 
estimation of the values at the steady concentration are attained by means of Eqs.([T7])- 



(18). Finally combining the relationships (19) with the definitions for {A, C, A', C'} in 



Eqs.©-®, we can then estimate these parameters as follows 

A = -2 (Y A -Y B )>0 

C ry QA,eros. i nB,evos. . rYj4,redist. . r(_B, redist. 

— ^\ c A,0^XY ~t~ c B,0JxY ~r c Afld X Y + c B,0 5 xy 

2 ' ■ (20) 

A' = -— (c Aib Y B + c Btb Y A ) > 



A 

2 



z^tl ( rr^,redist. riS,redist. 

° — ~~Z y C B,bC A flO xy ~~ c A,bCB,0^X,Y 



We see that fully four out of six of the BS parameters are accessible to MD studies. 
These are notable for a number of reasons. Of significant general interest, these happen to 
be precisely the parameters that govern the behavior of the longest wave perturbations to 



the surface height, via the finite-wavelength requirement in Eq.([T2]). Thus, the results of 
MD simulation may be used to predict whether an irradiated binary system has long-wave 
or finite-wavelength instabilities. Of more specific interest here, these parameters include 
both C and C, and so may be used to test specific requirements of the BS theory 
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IV. MOLECULAR-DYNAMICS ESTIMATION OF COEFFICIENTS 



A. Simulation Environment 

Using molecular dynamics simulations, we have estimated the quantities in Eq.([T6|) for 
the case GaSb irradiated by Ar ions at 250 eV. Amorphous GaSb was created by annealing a 
cubic piece of initially crystalline GaSb, consisting of 1.1 x 10 4 atoms. In these simulations 
periodic boundary conditions were used. Copies of the obtained amorphous GaSb were 
joined together (3 in the x- and y-directions, and 2 in the ^-direction) to form a target 
consisting of nearly 2.0 x 10 5 atoms. The surface of the target was finally relaxed at K, 
with periodicity only in the x- and y-directions. In the bombardment simulations periodic 
boundary conditions were used in the x- and y-directions. The simulation time was 250 ps, 
which is rather long, but turned out to be necessary for the post-impact atomic movements 
to seize. 

The inter-atomic potentials used for Ga-Sb interactions was that of Powell et al U34H . 
The (rarely occurring) pure Ga interactions were the same as those in the GaAs potential 
developed by us 113511 . To obtain the Sb potential, we note that As and Sb have the same 
crystal structure. Hence, the As potential parameter D e [1331] was rescaled with the differ- 
ence in cohesive energies of the two elements, and the parameter r with the difference in 
lattice constants. This gave a cohesive energy of 2.72 eV/atom and a density of 7.28 g/cm 3 
for pure Sb, in reasonable agreement with the experimental values of 2.72 eV/atom and 
6.74 g/cm 3 . We emphasize that since regions of pure Ga or Sb should not form in the 
simulated system, the results are not expected to be sensitive to small inaccuracies in the 
description of the pure elements. 

At each of nine different irradiation angles, we simulated 300 independent Ar ion im- 
pacts onto the target just constructed. Each ion hit the target at a random location, and for 
off normal incidence, arrived from a random azimuthal angle. As discussed in Ref. lT23l . 
two sources of noise were first filtered from the data before the moments were measured. 
First, despite the careful preparation of the target, some residual stresses remain, leading 
to displacements during the simulation that are due to the pre-existing state of the target, 
rather than to the impact. Second, the periodic nature of the cell allows small coherent 
shears to occur, on the scale of the entire cell. These are both filtered by constructing an 
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Figure 1. Moments and coefficients by component, (a) The decomposition of the zeroth moment 
M®, indicates that Ga is preferentially sputtered from a homogeneous amorphous film at 50/50 
initial concentrations, (b) The decomposition of the first erosive moment M e r2 s . indicates that 
Ga endures a stronger curvature-dependent sputtering, (c) The decomposition of the first redis- 
tributive moment M^\ ist indicates the preferential redistribution is not a dominant effect, (d) The 
erosive component of Sx for both species, (e) The redistributive component of Sx for both species, 
(f) A comparison of net erosive and redistributive contributions to the net coefficient Sx- (g)-(i) 
The corresponding plots for Sy (for completeness). As has been seen for pure silicon, the stabilizing 
redistributive contribution overwhelms the destabilizing erosive contribution at low angles. 



annulus around the impact cite, and calculating both the average lab-frame displacement 
of each atom, as well as the average shear in a co-ordinate system aligned with the incom- 
ing ion direction. Both averages were then subtracted from the displacement field, after 



which moments were calculated according to Eq.(|14|) 
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B. Results 



The results of our simulations are shown in FigjTJ We first calculate the moments ac- 
cording to Eq.( 14), and plot them along with three-term Fourier fits with expected symme- 
tries imposed (see Ref. ll23l0 . Our FigJT^ indicates that, for an initial fresh, unsegregated 



50/50 target, Ga is preferentially sputtered relative to Sb (see Section IV C below for a dis- 



cussion of this result); therefore we set A = Ga and B = Sb in the framework of Section III 



A comparison of Figure [TJd and [T]: supports our hypothesis that redistributive components 
of the first moment are of larger magnitude than the erosive components. In addition, the 
redistributive moments appear strikingly similar between species, indicating that preferen- 
tial redistribution is unlikely to be a significant effect for this system. To confirm that these 



impressions are not altered under the derivatives in Eqs.(|T6|), we also plot the coefficients 
found in that equation in Figures [TJd-i. We see, as in the case of pure silicon, an expected 
positive value for the net values of S x and S Y at normal incidence, when summed over all 
species and effects. Furthermore, we are able to provide the first prediction of the effect 
of preferential redistribution, and have found it unlikely to be important. 

Using the data shown in Figure [TJ and the definitions of the coefficients from Eqs.([T6]), 
we obtain the values 



S' 



S 



y^Ga 

ySb 

Ga,eros. 
X,Y 

Sb,eros. 



8=0 



S 



S 



X,Y 

Ga,redist. 
X,Y 

Sb,redist. 
X,Y 



8=0 



8=0 



8=0 



8=0 



-0.0172 x 7 
-0.0102 x 7 
-0.269 x L 



nm 

s 

nm 



nm 



-0.176 x I, 



nm 



(21) 



o 



1.11 x L 



nm 



o 



1.22 x I, 



nm 



o 



where I is again the flux in the beam direction, as in Eq.(|T6|). From Eq. (12) of Ref. HT4H ^ 
together with the relationships in Eq|l9] for Y Ga and Y Sb , we can obtain the predicted 
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steady concentrations of Ga and Sb to be 



CGa,0 = 0.37 

csb,o = 0.63 



Then, following the relationships (20), we convert these values into estimates for {A, C, A', C'} 
of 



A 
C 
A' 

a 



0.014 x I 
1.94 x J„ 
.0091 x J 
-.12 x I 



nm 



nm 



s 

nm 



where we have taken the value A « 3 nm from the irradiation of pure silicon at 250 eV 
11361 . These estimations exhibit three salient features: (1) the parameter C is positive, 
consistent with findings for pure materials, and contrary to the assumptions of BS theory; 
(2) the parameter group C'A in Eq.Q (the only place that C occurs) is much smaller than 
the similarly-dimensioned group A'C, indicating that preferential redistribution is likely 
not an important effect; and (3) the parameter group (A'C — C'A) is positive, indicating 
that long wave perturbations are stabilized. The first two findings are contrary to the 
assumptions of BS theory. Together with the third, and according to the discussion in 
Section |II C\ they imply that, in fact, that theory should predict stable, flat surfaces for 
irradiation of GaSb. 



C. Limitations. 

We want to acknowledge the following limitations in our approach. 



First, the definitions for the coefficients Y z and S z ' typc in Eq.(16), like the correspond- 
ing definitions in Ref. 11231 from which they are generalized, are the result of performing 
simulated impacts on an initially flat target surface. As such, they represent simplifications 
of the more general results in Ref. ||2~7ll for curved surfaces. Obtaining the latter would, 
unfortunately, require simulation of impacts not only over many angles, but also on many 
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different types of curved surface, which is outside our current computational capacity. It 
is our belief that this simplification does not significantly change the relative strength of 
erosion vs. redistribution. 

Second, our simulations were performed only at the 50/50 concentration for a fresh, 
homogeneous target, and values of coefficients at the steady concentration were inferred 
by means of the linear approximations in Eqs. (JT7]) - (f5~8|) . Like the flat-surface simplifi- 



cation, this avoids the simulation of targets at a variety of concentrations. Were more 
accurate estimates desired, such simulations could of course be performed. However, in 
addition to the computational cost of the simulations themselves, this would require the 
construction of many targets, which is in fact the most time-consuming part of our study. 
In addition, the nature of the irradiated target away from the 50/50 concentration is not 
currently well-known. 

Third, our results indicate that Ga is sputtered preferentially relative to Sb for a fresh, 
unsegregated target at the 50/50 concentration. This would imply the rapid enrichment 
of Sb near the surface, in direct contrast to experimental observations of excess Ga at the 
surface of this material, (more recent in situ observations made in an environment devoid 
of oxygen show that, while not as strong in the absence of oxygen, Ga enrichment is still 
observed ||29ll .) Because of the disparate experimental results, our MD result was checked 
against the results of SRIM, which predicted much more similar yields of the two species, 
but still with a 2% enhanced yield for Ga over Sb. Both the MD and the BCA thus give a 
result consistent with the more-efficient transfer of energy in Ar/Ga collisions relative to 
that in Ar/Sb collisions. It may thus be that actual experimental conditions differ from the 
ones we have simulated in ways that still need to be understood. 

Bringing the simulated results into agreement with experiment will require a better 
understanding of precisely what is happening during the early stages of irradiation, and 
the construction of an atomistic target reflecting that understanding. Nevertheless, an- 
ticipate that our principal conclusions are likely to remain valid for a more accurate set 
of targets. For instance, we found the redistributive component of the coefficient C - of 
central importance to the Bradley-Shipman theory - to be around five times greater than 
the erosive component. Because C represents the effects of net erosion and redistribution, 
it therefore should not be especially sensitive to the exact composition of the target, and 
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this dominance should be expected to continue. Similarly, although the groups AC and 
CA' contain the effects of preferential displacements, and may be expected to vary more 
with concentration, we find that CA' is more than ten times greater than AC - a difficult 
ratio to reverse under any circumstances. 

V. CONCLUSIONS 

Subject to the limitations just described, we have shown that four parameters of the 
Bradley-Shipman theory and other similar theories are accessible to Molecular Dynamics 
simulation by means of an extension, to binary materials, of the theory of crater functions 
described in 11231 |27ll . These parameters are A, C, A', and C, describing net sputtering 
and redistribution (C), preferential sputtering (A), preferential redistribution (C), and 
material replenishment (A') . Remarkably, these are exactly the four parameters needed 
to calculate the critical parameter group (A'C — CA), which governs the behavior of long 
waves. The methods demonstrated here are general, and should therefore be applicable 
to a wide variety of coupled PDE models for irradiated binary systems [|T3l[T4H , enabling 
the prediction of this important distinguishing characteristic among a variety of pattern- 
forming systems. 

Subsequently, we have performed the first estimate of these four parameters for the 
specific system of low-energy irradiation of GaSb by Ar, in which hexagonal dot arrays were 
first observed. Our main finding for this system is that the parameter C, associated with net 
erosion vs. net redistribution, is positive, indicating the dominance of redistribution over 
erosion, as was predicted by MD for pure Si [23], and shown experimentally for both pure 
Si and Ge |fT9l |20H . Additionally, we saw that the parameter 7 associated with preferential 
redistribution, hypothesized to be an important physical effect in ordered dot formation, is 
too small to play an important role, at least for this particular system. An interesting topic 
for future study will be to determine whether this tends to be true in general. 

Both the positive value of C, and the relative unimportance of C , are contrary to the 
finite-wavelength requirements of the Bradley-Shipman theory. In fact, for the parame- 
ter values we have estimated, that theory would predict entirely smooth surfaces for low- 
energy normal-incidence irradiation of GaSb. Especially if the dominance of the redistribu- 
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tive component of the collision cascade extends to a wide range of other two-component 
materials (in other words, if it is a generic property of the kinetically-dominated collision 
cascade), then this finding provides strong motivation for the consideration of alterna- 
tive physical mechanisms to explain the ordered patterns observed in these systems. We 
investigate such a mechanism elsewhere Il37ll . 
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